home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / cblas / source_hemm.h < prev    next >
Encoding:
C/C++ Source or Header  |  2001-05-14  |  7.4 KB  |  205 lines
  1. /* blas/source_hemm.h
  2.  * 
  3.  * Copyright (C) 2001 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. {
  21.   INDEX i, j, k;
  22.   INDEX n1, n2;
  23.   int uplo, side;
  24.  
  25.   const BASE alpha_real = CONST_REAL0(alpha);
  26.   const BASE alpha_imag = CONST_IMAG0(alpha);
  27.  
  28.   const BASE beta_real = CONST_REAL0(beta);
  29.   const BASE beta_imag = CONST_IMAG0(beta);
  30.  
  31.   if ((alpha_real == 0.0 && alpha_imag == 0.0)
  32.       && (beta_real == 1.0 && beta_imag == 0.0))
  33.     return;
  34.  
  35.   if (Order == CblasRowMajor) {
  36.     n1 = M;
  37.     n2 = N;
  38.     uplo = Uplo;
  39.     side = Side;
  40.   } else {
  41.     n1 = N;
  42.     n2 = M;
  43.     uplo = (Uplo == CblasUpper) ? CblasLower : CblasUpper;
  44.     side = (Side == CblasLeft) ? CblasRight : CblasLeft;
  45.   }
  46.  
  47.   /* form  y := beta*y */
  48.   if (beta_real == 0.0 && beta_imag == 0.0) {
  49.     for (i = 0; i < n1; i++) {
  50.       for (j = 0; j < n2; j++) {
  51.     REAL(C, ldc * i + j) = 0.0;
  52.     IMAG(C, ldc * i + j) = 0.0;
  53.       }
  54.     }
  55.   } else if (!(beta_real == 1.0 && beta_imag == 0.0)) {
  56.     for (i = 0; i < n1; i++) {
  57.       for (j = 0; j < n2; j++) {
  58.     const BASE Cij_real = REAL(C, ldc * i + j);
  59.     const BASE Cij_imag = IMAG(C, ldc * i + j);
  60.     REAL(C, ldc * i + j) = beta_real * Cij_real - beta_imag * Cij_imag;
  61.     IMAG(C, ldc * i + j) = beta_real * Cij_imag + beta_imag * Cij_real;
  62.       }
  63.     }
  64.   }
  65.  
  66.   if (alpha_real == 0.0 && alpha_imag == 0.0)
  67.     return;
  68.  
  69.   if (side == CblasLeft && uplo == CblasUpper) {
  70.  
  71.     /* form  C := alpha*A*B + C */
  72.  
  73.     for (i = 0; i < n1; i++) {
  74.       for (j = 0; j < n2; j++) {
  75.     const BASE Bij_real = CONST_REAL(B, ldb * i + j);
  76.     const BASE Bij_imag = CONST_IMAG(B, ldb * i + j);
  77.     const BASE temp1_real = alpha_real * Bij_real - alpha_imag * Bij_imag;
  78.     const BASE temp1_imag = alpha_real * Bij_imag + alpha_imag * Bij_real;
  79.     BASE temp2_real = 0.0;
  80.     BASE temp2_imag = 0.0;
  81.     {
  82.       const BASE Aii_real = CONST_REAL(A, i * lda + i);
  83.       /* const BASE Aii_imag = 0.0; */
  84.       REAL(C, i * ldc + j) += temp1_real * Aii_real;
  85.       IMAG(C, i * ldc + j) += temp1_imag * Aii_real;
  86.     }
  87.     for (k = i + 1; k < n1; k++) {
  88.       const BASE Aik_real = CONST_REAL(A, i * lda + k);
  89.       const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
  90.       const BASE Bkj_real = CONST_REAL(B, ldb * k + j);
  91.       const BASE Bkj_imag = CONST_IMAG(B, ldb * k + j);
  92.       REAL(C, k * ldc + j) += Aik_real * temp1_real - (-Aik_imag) * temp1_imag;
  93.       IMAG(C, k * ldc + j) += Aik_real * temp1_imag + (-Aik_imag) * temp1_real;
  94.       temp2_real += Aik_real * Bkj_real - Aik_imag * Bkj_imag;
  95.       temp2_imag += Aik_real * Bkj_imag + Aik_imag * Bkj_real;
  96.     }
  97.     REAL(C, i * ldc + j) += alpha_real * temp2_real - alpha_imag * temp2_imag;
  98.     IMAG(C, i * ldc + j) += alpha_real * temp2_imag + alpha_imag * temp2_real;
  99.       }
  100.     }
  101.  
  102.   } else if (side == CblasLeft && uplo == CblasLower) {
  103.  
  104.     /* form  C := alpha*A*B + C */
  105.  
  106.     for (i = 0; i < n1; i++) {
  107.       for (j = 0; j < n2; j++) {
  108.     const BASE Bij_real = CONST_REAL(B, ldb * i + j);
  109.     const BASE Bij_imag = CONST_IMAG(B, ldb * i + j);
  110.     const BASE temp1_real = alpha_real * Bij_real - alpha_imag * Bij_imag;
  111.     const BASE temp1_imag = alpha_real * Bij_imag + alpha_imag * Bij_real;
  112.     BASE temp2_real = 0.0;
  113.     BASE temp2_imag = 0.0;
  114.     for (k = 0; k < i; k++) {
  115.       const BASE Aik_real = CONST_REAL(A, i * lda + k);
  116.       const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
  117.       const BASE Bkj_real = CONST_REAL(B, ldb * k + j);
  118.       const BASE Bkj_imag = CONST_IMAG(B, ldb * k + j);
  119.       REAL(C, k * ldc + j) += Aik_real * temp1_real - (-Aik_imag) * temp1_imag;
  120.       IMAG(C, k * ldc + j) += Aik_real * temp1_imag + (-Aik_imag) * temp1_real;
  121.       temp2_real += Aik_real * Bkj_real - Aik_imag * Bkj_imag;
  122.       temp2_imag += Aik_real * Bkj_imag + Aik_imag * Bkj_real;
  123.     }
  124.     {
  125.       const BASE Aii_real = CONST_REAL(A, i * lda + i);
  126.       /* const BASE Aii_imag = 0.0; */
  127.       REAL(C, i * ldc + j) += temp1_real * Aii_real;
  128.       IMAG(C, i * ldc + j) += temp1_imag * Aii_real;
  129.     }
  130.     REAL(C, i * ldc + j) += alpha_real * temp2_real - alpha_imag * temp2_imag;
  131.     IMAG(C, i * ldc + j) += alpha_real * temp2_imag + alpha_imag * temp2_real;
  132.       }
  133.     }
  134.  
  135.   } else if (side == CblasRight && uplo == CblasUpper) {
  136.  
  137.     /* form  C := alpha*B*A + C */
  138.  
  139.     for (i = 0; i < n1; i++) {
  140.       for (j = 0; j < n2; j++) {
  141.     const BASE Bij_real = CONST_REAL(B, ldb * i + j);
  142.     const BASE Bij_imag = CONST_IMAG(B, ldb * i + j);
  143.     const BASE temp1_real = alpha_real * Bij_real - alpha_imag * Bij_imag;
  144.     const BASE temp1_imag = alpha_real * Bij_imag + alpha_imag * Bij_real;
  145.     BASE temp2_real = 0.0;
  146.     BASE temp2_imag = 0.0;
  147.     {
  148.       const BASE Ajj_real = CONST_REAL(A, j * lda + j);
  149.       /* const BASE Ajj_imag = 0.0; */
  150.       REAL(C, i * ldc + j) += temp1_real * Ajj_real;
  151.       IMAG(C, i * ldc + j) += temp1_imag * Ajj_real;
  152.     }
  153.     for (k = j + 1; k < n2; k++) {
  154.       const BASE Ajk_real = CONST_REAL(A, j * lda + k);
  155.       const BASE Ajk_imag = CONST_IMAG(A, j * lda + k);
  156.       const BASE Bik_real = CONST_REAL(B, ldb * i + k);
  157.       const BASE Bik_imag = CONST_IMAG(B, ldb * i + k);
  158.       REAL(C, i * ldc + k) += temp1_real * Ajk_real - temp1_imag * Ajk_imag;
  159.       IMAG(C, i * ldc + k) += temp1_real * Ajk_imag + temp1_imag * Ajk_real;
  160.       temp2_real += Bik_real * Ajk_real - Bik_imag * (-Ajk_imag);
  161.       temp2_imag += Bik_real * (-Ajk_imag) + Bik_imag * Ajk_real;
  162.     }
  163.     REAL(C, i * ldc + j) += alpha_real * temp2_real - alpha_imag * temp2_imag;
  164.     IMAG(C, i * ldc + j) += alpha_real * temp2_imag + alpha_imag * temp2_real;
  165.       }
  166.     }
  167.  
  168.   } else if (side == CblasRight && uplo == CblasLower) {
  169.  
  170.     /* form  C := alpha*B*A + C */
  171.  
  172.     for (i = 0; i < n1; i++) {
  173.       for (j = 0; j < n2; j++) {
  174.     const BASE Bij_real = CONST_REAL(B, ldb * i + j);
  175.     const BASE Bij_imag = CONST_IMAG(B, ldb * i + j);
  176.     const BASE temp1_real = alpha_real * Bij_real - alpha_imag * Bij_imag;
  177.     const BASE temp1_imag = alpha_real * Bij_imag + alpha_imag * Bij_real;
  178.     BASE temp2_real = 0.0;
  179.     BASE temp2_imag = 0.0;
  180.     for (k = 0; k < j; k++) {
  181.       const BASE Ajk_real = CONST_REAL(A, j * lda + k);
  182.       const BASE Ajk_imag = CONST_IMAG(A, j * lda + k);
  183.       const BASE Bik_real = CONST_REAL(B, ldb * i + k);
  184.       const BASE Bik_imag = CONST_IMAG(B, ldb * i + k);
  185.       REAL(C, i * ldc + k) += temp1_real * Ajk_real - temp1_imag * Ajk_imag;
  186.       IMAG(C, i * ldc + k) += temp1_real * Ajk_imag + temp1_imag * Ajk_real;
  187.       temp2_real += Bik_real * Ajk_real - Bik_imag * (-Ajk_imag);
  188.       temp2_imag += Bik_real * (-Ajk_imag) + Bik_imag * Ajk_real;
  189.     }
  190.     {
  191.       const BASE Ajj_real = CONST_REAL(A, j * lda + j);
  192.       /* const BASE Ajj_imag = 0.0; */
  193.       REAL(C, i * ldc + j) += temp1_real * Ajj_real;
  194.       IMAG(C, i * ldc + j) += temp1_imag * Ajj_real;
  195.     }
  196.     REAL(C, i * ldc + j) += alpha_real * temp2_real - alpha_imag * temp2_imag;
  197.     IMAG(C, i * ldc + j) += alpha_real * temp2_imag + alpha_imag * temp2_real;
  198.       }
  199.     }
  200.  
  201.   } else {
  202.     BLAS_ERROR("unrecognized operation");
  203.   }
  204. }
  205.